package scu.maqiang.fes;

import static scu.maqiang.numeric.Constants.GP3;
import static scu.maqiang.numeric.Constants.GW3;
import static scu.maqiang.numeric.MVO.matmul;
import static scu.maqiang.numeric.MVO.mul;
import static scu.maqiang.numeric.MVO.transpose;

import java.util.Arrays;

import scu.maqiang.mesh.MatrixFunc;
import scu.maqiang.mesh.Mesh2Q8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.CG;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class FES2Q82 extends FES {
	private static final double[] GP4 = null;
	private static final int[] GW4 = null;

	public FES2Q82(Mesh2Q8 mesh2D) {
		super(mesh2D, 2);
		eleStiff = FES2Q82::elementElasticityStiff;
		eleMass = FES2Q82::elementElasticityMass;
		eleSource = FES2Q82::elementVolumeForce;
		eleFlux = FES2T32::elemenSurfaceForce;
	}
	
	/**
	 * ���㵯�����ⵥԪ�նȾ���
	 * 
	 * @param xy
	 *            3x2ά��Ԫ�������ÿ�б�ʾ����x��y����
	 * @param ENu
	 *            2x3ά���Գ���
	 * @param tp
	 * @param K
	 */
	public static void elementElasticityStiff(int[] ele, double[][] xy, double[][] ENu, BVPType tp, double[][] K) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double xi, et, detJ = 0.0, wxarea = 0.0, c;
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] v = new double[8];
		double[] E = ENu[0];
		double[] Nu = ENu[1];
		double mE, mNu;
		double[][] KK = new double[16][16];
		
		double[][] B = null;
		MVO.fill(K, 0.0);
		switch (tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
			B = new double[3][16];
			for(int i = 0; i < 3; i++) {
				xi = GP3[i];
				for(int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					mE = MVO.dot_product(v, E);
					mNu = MVO.dot_product(v, Nu);
					double[][] D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
					wxarea = GW3[i] * GW3[j] * detJ;
					for(int ib = 0; ib < 8; ib++) {
						B[0][2 * ib] = B[2][2 * ib + 1] = dx[ib];
						B[1][2 * ib + 1] = B[2][2 * ib] = dy[ib];
					}
					mul(matmul(transpose(B), matmul(D, B)), wxarea, KK);
					for(int ii = 0; ii < 16; ii++) {
						for(int jj = 0; jj < 16; jj++) {
							K[ii][jj] += KK[ii][jj];
						}
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_COMMON:
			B = new double[3][16];
			for(int i = 0; i < 3; i++) {
				xi = GP3[i];
				for(int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					wxarea = GW3[i] * GW3[j] * detJ;
					for(int ib = 0; ib < 8; ib++) {
						B[0][2 * ib] = B[2][2 * ib + 1] = dx[ib];
						B[1][2 * ib + 1] = B[2][2 * ib] = dy[ib];
					}
					mul(matmul(transpose(B), matmul(ENu, B)), wxarea, KK);
					for(int ii = 0; ii < 16; ii++) {
						for(int jj = 0; jj < 16; jj++) {
							K[ii][jj] += KK[ii][jj];
						}
					}
				}
			}
			break;
		case AXISYMMETRIC:
			B = new double[4][16];
			for(int i = 0; i < 3; i++) {
				xi = GP3[i];
				for(int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					mE = MVO.dot_product(v, E);
					mNu = MVO.dot_product(v, Nu);
					double[][] D = Materials.ElasticConstituteMatrix(mE, mNu, tp);
					double r0 = MVO.dot_product(x, v);
					wxarea = GW3[i] * GW3[j] * detJ * r0;
					for(int ib = 0; ib < 8; ib++) {
						B[0][2 * ib] = B[3][2 * ib + 1] = dx[ib];
						B[1][2 * ib + 1] = B[3][2 * ib] = dy[ib];
						B[2][2 * ib] = v[ib] / r0;
					}
					mul(matmul(transpose(B), matmul(D, B)), wxarea, KK);
					for(int ii = 0; ii < 16; ii++) {
						for(int jj = 0; jj < 16; jj++) {
							K[ii][jj] += KK[ii][jj];
						}
					}
				}
			}
			break;
		case CONSTITUTE_MATRIX_AXIS:
			B = new double[4][16];
			for(int i = 0; i < 3; i++) {
				xi = GP3[i];
				for(int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					double r0 = MVO.dot_product(x, v);
					wxarea = GW3[i] * GW3[j] * detJ * r0;
					for(int ib = 0; ib < 8; ib++) {
						B[0][2 * ib] = B[3][2 * ib + 1] = dx[ib];
						B[1][2 * ib + 1] = B[3][2 * ib] = dy[ib];
						B[2][2 * ib] = v[ib] / r0;
					}
					mul(matmul(transpose(B), matmul(ENu, B)), wxarea, KK);
					for(int ii = 0; ii < 16; ii++) {
						for(int jj = 0; jj < 16; jj++) {
							K[ii][jj] += KK[ii][jj];
						}
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Wrong BVPType!");
		}
//		System.out.println(MVO.toString(K));
//		System.exit(0);
	}

	public static void elementElasticityMass(int[] ele, double[][] xy, double[][] rho, BVPType tp, double[][] eleMatrix) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] v = new double[8];
		double xi, et, detJ, wxarea, c;
		MVO.fill(eleMatrix, 0.0);
		switch (tp) {
		case CONSIST_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(rho[0], v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						for (int sj = 0; sj < 8; sj++) {
							eleMatrix[2 * si][2 * sj] += v[si] * v[sj] * wxarea;
							eleMatrix[2 * si + 1][2 * sj + 1] += v[si] * v[sj] * wxarea;
						}
					}
				}
			}
			break;
		case LUMP_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(rho[0], v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						eleMatrix[2 * si][2 * si] += v[si] * wxarea;
						eleMatrix[2 * si + 1][2 * si + 1] += v[si] * wxarea;
					}
				}
			}
			break;
		case AXIS_CONSIST_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(rho[0], v) * MVO.dot_product(x, v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						for (int sj = 0; sj < 8; sj++) {
							eleMatrix[2 * si][2 * sj] += v[si] * v[sj] * wxarea;
							eleMatrix[2 * si + 1][2 * sj + 1] += v[si] * v[sj] * wxarea;
						}
					}
				}
			}
			break;
		case AXIS_LUMP_MASS:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					c = MVO.dot_product(rho[0], v) * MVO.dot_product(x, v);
					wxarea = GW3[i] * GW3[j] * detJ * c;
					for (int si = 0; si < 8; si++) {
						eleMatrix[2 * si][2 * si] += v[si] * wxarea;
						eleMatrix[2 * si + 1][2 * si + 1] += v[si] * wxarea;
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Illegal BVPType!");
		}
	}

	public static void elementVolumeForce(int[] ele, double[][] xy, double[][] f, BVPType tp, double[] eleF) {
		double[] x = MVO.col(xy, 0);
		double[] y = MVO.col(xy, 1);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] v = new double[8];
		double xi, et, detJ, wxarea, fu, fv;
		Arrays.fill(eleF, 0.0);
		switch (tp) {
		case COMMON:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					fu = MVO.dot_product(v, f[0]);
					fv = MVO.dot_product(v, f[1]);
					wxarea = GW3[i] * GW3[j] * detJ;
					for (int si = 0; si < 8; si++) {
						ele[2 * si] += v[si]  * fu *wxarea;
						ele[2 * si + 1] += v[si] * fv * wxarea;
					}
				}
			}
			break;
		case AXISYMMETRIC:
			for (int i = 0; i < 3; i++) {
				xi = GP3[i];
				for (int j = 0; j < 3; j++) {
					et = GP3[j];
					detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
					fu = MVO.dot_product(f[0], v) ;
					fv = MVO.dot_product(v, f[1]);
					wxarea = GW3[i] * GW3[j] * detJ * MVO.dot_product(x, v);
					for (int si = 0; si < 8; si++) {
						ele[2 * si] += v[si] * fu * wxarea;
						ele[2 * si + 1] += v[si] * fv * wxarea;
					}
				}
			}
			break;
		default:
			throw new IllegalArgumentException("Illegal BVPType!");
		}
	}

	public void assemblePressure(ScalarFunc[] coefFunc, double[] param, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(coefFunc, param, FES2T62::elementPressure, tp, RHS, bdLabel);
	}

	public void assemblePressure(double[] constCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(constCoef, FES2T62::elementPressure, tp, RHS, bdLabel);
	}
	
	public void assemblePressure(double[][] vecCoef, BVPType tp, double[] RHS, int... bdLabel) {
		assembleBoundaryVector(vecCoef, FES2T62::elementPressure, tp, RHS, bdLabel);
	}

	public void assembleHeatLoad(ScalarFunc[] coef, double[] param, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, param, vecCoef, FES2Q82::elementHeatLoad, tp, RHS);
	}

	public void assembleHeatLoad(double[] coef, double[][] vecCoef, BVPType tp, double[] RHS) {
		assembleGlobalVector(coef, vecCoef, FES2Q82::elementHeatLoad, tp, RHS);
	}

	public static void elementHeatLoad(int[] ele, double[][] coord, double[][] coef, BVPType tp, double[] eleLoad) {
		double[] x = MVO.col(coord, 0);
		double[] y = MVO.col(coord, 1);
		double[] dx = new double[8];
		double[] dy = new double[8];
		double[] v = new double[8];
		double xi, et, detJ, wxarea;
		Arrays.fill(eleLoad, 0.0);
	    double mx;
	    switch(tp) {
		case PLANESTRAIN:
		case PLANESTRESS:
		case COMMON:
	    	for(int i = 0; i < 3; i++) {
	    		xi = GP3[i];
	    		for(int j = 0; j < 3; j++) {
	    			et = GP3[j];
	    			detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
	    			wxarea = GW3[i] * GW3[j] * detJ * MVO.dot_product(v, coef[0]) * MVO.dot_product(v, coef[1]);
	    			for(int ii = 0; ii < 8; ii++) {
	    				eleLoad[2 * ii] += dx[ii] * wxarea;
	    				eleLoad[2 * ii + 1] += dy[ii] * wxarea;
	    			}
	    		}
	    	}
	    	break;
	    case CONSTITUTE_MATRIX_COMMON:
	    	break;
		case AXISYMMETRIC:
	    	for(int i = 0; i < 4; i++) {
	    		xi = GP4[i];
	    		for(int j = 0; j < 4; j++) {
	    			et = GP4[j];
	    			detJ = CG.ISOMap2D(CG::Q8ShapeFunction, x, y, xi, et, v, dx, dy);
	    			wxarea = GW4[i] * GW4[j] * detJ * MVO.dot_product(v, coef[0]) * MVO.dot_product(v, coef[1]);
	                mx = MVO.dot_product(x, v);
	    			for(int ii = 0; ii < 8; ii++) {
	    				eleLoad[2 * ii] += (dx[ii]  * mx + v[ii]) * wxarea;
	    				eleLoad[2 * ii + 1] += dy[ii] * mx * wxarea;
	    			}
	    		}
	    	}
	    	break;
		case CONSTITUTE_MATRIX_AXIS:
			break;
	    default:
	    	throw new IllegalArgumentException("Illegal BVPType!");	
	    }
	}

	public void computeNodalStrain(double[] sol, double[][] strain) {
//		int len = sol.length;
//		ParamCheck.checkEqual(len, this.ndof);
//		int n = len / 2;
//		ParamCheck.checkEqual(n, strain[0].length);
//		ParamCheck.checkEqual(n, strain[1].length);
//		ParamCheck.checkEqual(n, strain[2].length);
//		double[] angle = new double[3];
//		double[][] coord = new double[3][2];
//		double[][] grad = new double[3][2];
//		double[] angleCount = new double[this.GetNdof()];
//		double dxU, dyU, dxV, dyV;
//		double[] eleU = new double[3];
//		double[] eleV = new double[3];
//		for (int i = 0; i < mesh.getNt(); i++) {
//			int[] ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			CG.gradient2DT3(coord, grad);
//			CG.angle2DT3(coord, angle);
//			for (int j = 0; j < 3; j++) {
//				eleU[j] = sol[2 * ele[j]];
//				eleV[j] = sol[2 * ele[j] + 1];
//			}
//			double[] dx = MVO.col(grad, 0);
//			double[] dy = MVO.col(grad, 1);
//			dxU = MVO.dot_product(eleU, dx);
//			dyU = MVO.dot_product(eleU, dy);
//			dxV = MVO.dot_product(eleV, dx);
//			dyV = MVO.dot_product(eleV, dy);
//			for (int j = 0; j < 3; j++) {
//				strain[0][ele[j]] += dxU * angle[j];
//				strain[1][ele[j]] += dyV * angle[j];
//				strain[2][ele[j]] += (dxV + dyU) * angle[j];
//				angleCount[ele[j]] += angle[j];
//			}
//		}
//		for (int i = 0; i < n; i++) {
//			strain[0][i] /= angleCount[i];
//			strain[1][i] /= angleCount[i];
//			strain[2][i] /= angleCount[i];
//		}
	}

	public void computeElementStrain(double[] sol, double[][] strain) {
//		int len = sol.length;
//		ParamCheck.checkEqual(len, this.ndof);
//		int n = len / 2;
//		int nt = mesh.getNt();
//		ParamCheck.checkEqual(nt, strain[0].length);
//		ParamCheck.checkEqual(nt, strain[1].length);
//		ParamCheck.checkEqual(nt, strain[2].length);
//		double[] angle = new double[3];
//		double[][] coord = new double[3][2];
//		double[][] grad = new double[3][2];
//		double dxU, dyU, dxV, dyV;
//		double[] eleU = new double[3];
//		double[] eleV = new double[3];
//		for (int i = 0; i < nt; i++) {
//			int[] ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			CG.gradient2DT3(coord, grad);
//			CG.angle2DT3(coord, angle);
//			for (int j = 0; j < 3; j++) {
//				eleU[j] = sol[2 * ele[j]];
//				eleV[j] = sol[2 * ele[j] + 1];
//			}
//			double[] dx = MVO.col(grad, 0);
//			double[] dy = MVO.col(grad, 1);
//			dxU = MVO.dot_product(eleU, dx);
//			dyU = MVO.dot_product(eleU, dy);
//			dxV = MVO.dot_product(eleV, dx);
//			dyV = MVO.dot_product(eleV, dy);
//			strain[0][i] = dxU;
//			strain[1][i] = dyV;
//			strain[2][i] = dxV + dyU;
//		}
	}

	public void computeNodalStress(double[] sol, double[][] stress) {

	}

	public void computeElementStress(double[] sol, double[][] stress) {

	}

	public void computeNodalStrainStressVonMisesStress(BVPType tp, double[] coef, double[] sol, double[][] strain,
			double[][] stress, double[] misesStress) {
//		int len = sol.length;
//		ParamCheck.checkEqual(len, this.ndof);
//		int n = len / 2;
//		ParamCheck.checkEqual(n, strain[0].length);
//		ParamCheck.checkEqual(n, strain[1].length);
//		ParamCheck.checkEqual(n, strain[2].length);
//		double[] angle = new double[3];
//		double[][] coord = new double[3][2];
//		double[][] grad = new double[3][2];
//		double[] angleCount = new double[this.GetNdof()];
//		double dxU, dyU, dxV, dyV;
//		double[] eleU = new double[3];
//		double[] eleV = new double[3];
//		double E = coef[0];
//		double Nu = coef[1];
//		double mu = E / (2 * (1 + Nu));
//		double lambda;
//		double ex, ey, exy, ez, tx, ty, tz, txy, ttheta;
//		double t1, t2, t3;
//		double meanStress;
//		if (tp == BVPType.PLANESTRESS) {
//			lambda = E * Nu / (1 - Nu * Nu);
//		} else {
//			lambda = E * Nu / ((1 - 2 * Nu) * (1 + Nu));
//		}
//		for (int i = 0; i < mesh.getNt(); i++) {
//			int[] ele = mesh.getElement(i);
//			mesh.getCoordPerElement(ele, coord);
//			CG.gradient2DT3(coord, grad);
//			CG.angle2DT3(coord, angle);
//			for (int j = 0; j < 3; j++) {
//				eleU[j] = sol[2 * ele[j]];
//				eleV[j] = sol[2 * ele[j] + 1];
//			}
//			double[] x = MVO.col(coord, 0);
//			double[] dx = MVO.col(grad, 0);
//			double[] dy = MVO.col(grad, 1);
//			dxU = MVO.dot_product(eleU, dx);
//			dyU = MVO.dot_product(eleU, dy);
//			dxV = MVO.dot_product(eleV, dx);
//			dyV = MVO.dot_product(eleV, dy);
//			// for(int j = 0; j < 3; j++) {
//			// strain[0][ele.idx[j]] += dxU * angle[j];
//			// strain[1][ele.idx[j]] += dyV * angle[j];
//			// strain[2][ele.idx[j]] += (dxV + dyU) * angle[j];
//			// if (BVPType.AXISYMMETRIC == tp) {
//			// strain[3][ele.idx[j]] += eleU[j] / x[j] * angle[j];
//			// stress[0][ele.idx[j]] += (lambda + 2 * mu) * strain[0][ele.idx[j]] + lambda *
//			// (strain[1][ele.idx[j]] + strain[3][ele.idx[j]]);
//			// stress[1][ele.idx[j]] += lambda * (strain[0][ele.idx[j]] +
//			// strain[3][ele.idx[j]]) + (lambda + 2 * mu) * strain[1][ele.idx[j]];
//			// stress[3][ele.idx[j]] += (lambda + 2 * mu) * strain[3][ele.idx[j]] + lambda *
//			// (strain[0][ele.idx[j]] + strain[1][ele.idx[j]]);
//			// } else {
//			// tx = (lambda + 2 * mu) * strain[0][ele.idx[j]] + lambda *
//			// strain[1][ele.idx[j]];
//			// ty = lambda * strain[0][ele.idx[j]] + (lambda + 2 * mu) *
//			// strain[1][ele.idx[j]];
//			// stress[0][ele.idx[j]] += tx;
//			// stress[1][ele.idx[j]] += ty;
//			// if (BVPType.PLANESTRAIN == tp) {
//			// stress[3][ele.idx[j]] += Nu * (tx + ty);
//			// } else if (BVPType.PLANESTRESS == tp) {
//			// strain[3][ele.idx[j]] += -Nu * (tx + ty) / E;
//			// }
//			// }
//			// stress[2][ele.idx[j]] += mu * strain[2][ele.idx[j]];
//			// angleCount[ele.idx[j]] += angle[j];
//			// }
//			// }
//			// for(int i = 0; i < n; i++) {
//			// strain[0][i] /= angleCount[i];
//			// strain[1][i] /= angleCount[i];
//			// strain[2][i] /= angleCount[i];
//			// strain[3][i] /= angleCount[i];
//			// stress[0][i] /= angleCount[i];
//			// stress[1][i] /= angleCount[i];
//			// stress[2][i] /= angleCount[i];
//			// stress[3][i] /= angleCount[i];
//			// }
//
//			for (int j = 0; j < 3; j++) {
//				strain[0][ele[j]] += dxU * angle[j];
//				strain[1][ele[j]] += dyV * angle[j];
//				strain[2][ele[j]] += (dxV + dyU) * angle[j];
//				if (BVPType.AXISYMMETRIC == tp) {
//					strain[3][ele[j]] += eleU[j] / x[j] * angle[j];
//				}
//				angleCount[ele[j]] += angle[j];
//			}
//		}
//		for (int i = 0; i < n; i++) {
//			strain[0][i] /= angleCount[i];
//			strain[1][i] /= angleCount[i];
//			strain[2][i] /= angleCount[i];
//			strain[3][i] /= angleCount[i];
//		}
//
//		for (int i = 0; i < n; i++) {
//			if (BVPType.AXISYMMETRIC == tp) {
//				stress[0][i] = (lambda + 2 * mu) * strain[0][i] + lambda * (strain[1][i] + strain[3][i]);
//				stress[1][i] = lambda * (strain[0][i] + strain[3][i]) + (lambda + 2 * mu) * strain[1][i];
//				stress[3][i] = (lambda + 2 * mu) * strain[3][i] + lambda * (strain[0][i] + strain[1][i]);
//			} else {
//				stress[0][i] = tx = (lambda + 2 * mu) * strain[0][i] + lambda * strain[1][i];
//				stress[1][i] = ty = lambda * strain[0][i] + (lambda + 2 * mu) * strain[1][i];
//				if (BVPType.PLANESTRAIN == tp) {
//					stress[3][i] = Nu * (tx + ty);
//				} else if (BVPType.PLANESTRESS == tp) {
//					strain[3][i] = -Nu * (tx + ty) / E;
//				}
//			}
//			stress[2][i] = mu * strain[2][i];
//			meanStress = (stress[0][i] + stress[1][i] + stress[3][i]) / 3;
//			t1 = stress[0][i] - meanStress;
//			t2 = stress[1][i] - meanStress;
//			t3 = stress[3][i] - meanStress;
//			misesStress[i] = 0.5 * (t1 * t1 + t2 * t2 + t3 * t3) + stress[2][i] * stress[2][i];
//			misesStress[i] = Math.sqrt(3.0 * misesStress[i]);
//		}
	}
}
